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o 

^vj Abstract. The purpose of this paper is to design an algorithm for the computation of the 

counterparty risk which is competitive in regards of a brute force "Monte-Carlo of Monte-Carlo" 

2 method (with nested simulations) . This is achieved using marked branching diffusions describing 

^^ a Galton-Watson random tree. Such an algorithm leads at the same time to a computation of 

^^ . the (bilateral) counterparty risk when we use the default-risky or counterparty-risklcss option 

values as mark-to-mar kct. Our method is illustrated by various numerical examples. 



Cr~' 1. Introduction 

^*^ The recent financial crisis has highhghted the importance of credit valuation adjustment when 

rt pricing derivative contracts. Bilateral counterparty risk is the risk that the issuer of a derivative 

C.|_H contract or the counterparty may default prior to the expiry and fail to make future payments. This 

market imperfection leads naturally for Markovian models to non-linear second-order parabolic 

partial differential equations (PDEs). More precisely, the non-linearity in the pricing equation 

affects none of the differential terms and depends on the positive part of the mark-to-markct value 

of the derivative upon default. We have a so-called semi-linear PDE. The numerical solution of 

^^ this equation is a formidable task that has attracted little attention from practitioners. For niulti- 

\^ asset portfolios, these PDEs which suffer from the curse of dimensionality cannot be solved with 

CO finite-difference schemes. We must rely on probabilistic methods. Up to now, it seems that a brute 

^N force intensive "Monte-Carlo of Monte-Carlo" method (with nested simulations) is the only tool 

ff") available for this task. 

O . 

fvi In this paper, we rely on new advanced non-linear Monte-Carlo methods for solving these semi- 

,— I linear PDEs. A first approach is to use the so-called first-order backward stochastic differential 

ILJ equations. Unfortunately, in practise this method requires the computation of conditional ex- 

pectations using regressions. Finding good quality regressors is notably difficult, especially for 
multi-asset portfolios. This leads us to introduce a new method based on branching diffusions de- 
scribing a marked Galton-Watson random tree. A similar algorithm can also be applied to obtain 
stochastic representations for solutions of a large class of semi-linear parabolic PDEs in which the 
non-linearity can be approximated by a polynomial function. 



> 



X 



2. Credit Valuation Adjustment 

2.1. Semi-linear PDEs. For completeness, we derive the PDE arising in counterparty risk valua- 
tion of a European derivative with a payoff i/i at maturity T. In short, depending on the (modeling) 
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choice of the mark-to-market value of the derivative upon defauh, we will get two types of semi- 
linear PDEs that can be schematically written as 

(1) dfU + Cu + vqu + riu'^ ~ 0, u{T,x)=ip{x) 
and 

(2) dtu + Cu + rau + riM + r2M+ = 0, w(T, x) = Tp{x) 

dtM + CM + r^M = 0, M{T,x) ^ ipix) 

C is the Ito generator of a multi-dimensional diffusion process and tj are arbitrary functions of t 
and X. 

2.2. PDE derivation. We assume the issuer is allowed to dynamically trade d underlying assets 
X. S M!J_. Additionally, in order to hedge his credit risk on the counterparty name, he can trade 
a default risky bond, denoted Pi . Furthermore, the values of the underlyings are not altered by 
the counterparty default which is modeled by a Poisson jump process. For the sake of simplicity, 
we consider a constant intensity. This assumption can be easily relaxed, in particular the intensity 
can follow an Ito diffusion. For use below, expressions with a subscript 2 denote counterparty 
quantities. We consider the case of a long position in a single derivative whose value we denote 
u. In practice netting agreements apply to the global mark-to-market value of a pool of derivative 
positions - u would then denote the aggregate value of these derivatives. The processes Xt,Pi 
satisfy under the risk-neutral measure P (we assume the market model is complete) 



dXt 

d2 



dP? 



P? 



rdt + a{t,Xt).dWt 

{r + X2)dt-dJf 

with Wt a d-dimensional Brownian motion, Jf a jump Poisson process with intensity A2 and r the 
interest rate. The no-arbitrage condition and the completeness of the market give that e~^*u{t, Xt) 
is a P-martingale, characterized by 

dtu + Cu + A2 (m — u) — ru = 

where C denotes the Ito generator of X and u the derivative value after the counterparty has 
defaulted. At the default event, u is given h^ 

u = RM+ - M- 

with AI the mark-to-market value of the derivative to be used in the unwinding of the position 
upon default and R the recovery rate. There is an ambiguity in the market about the convention 
for the mark-to-market value to be settled at default. There are two natural conventions (see [4] 
for discussions about the relevance of these conventions): The mark-to- market of the derivative is 
evaluated at the time of default with provision for counterparty risk or without. 

1. Provision for counterparty risk, M — u: 

(3) dtU + Cu-{1- R)X2U+ - ru = 0, u{T, x) = 7/;(a;) 

In the particular case when the payoff ^(a;) is negative, the solution is given by e'~^^'^~*'>Et.x['4'iXT)]- 

2. No provision for counterparty risk: 

(4) dtU + CU + X2 {RM+ - M- - u) - ru = 0, u{T, x) = i}{x) 
dtM + CM - rM = 0, M{T, x) = tl){x) 

^X = X+ -X- . 
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In the case of collateralized positions, counterparty risk applies to the variation of the mark-to- 
market value of the corresponding positions experienced over the time it takes to qualify a failure 
to pay margin as a default event - typically a few days. In the latter case, the non-linearity uf 
should be substituted with {ut — Ui+a)^ where A is this delay. We will come back to this situation 
in the last section (see remark [5^. 

By proper discounting and replacing u by —u for the sake of the presentation, these two PDEs can 
be cast into normal forms 

(5) dtu + Cu + P{u+ -u) :^ 0, u{T, x) = i:{x) : PDE2 

(6) dtu + Cu+^-^{{l~R)Et,,[^]+ + REt,4'4>]-u)^0, u{T, x) = ^{x) -.PDEl 

with /3 = A2(l — i?) G M+. It is interesting to note that a similar semi- linear PDE type ^ appears 
also in the pricing of American options. 

2.3. American options. The replication price of an American option with exercise payoff ip{x) 
satisfies a variational PDE: 

max (dtu + Cu, ip{x) — u) = 0, u{T, x) = ip{x) 

This PDE can be converted into a semi-linear PDE (see [2] for details): 

9tu + £u = 1^(^)>„£V(2;), u{T,x) ^ip{x) 

Stochastic representations of this equation lead to well-known early exercise premium formulas of 
American options. Our algorithm can also be applied to this non-linear PDE. It does not require 
regressions as in the well-known Longstaff-Schwartz method [10] or a "Monte-Carlo of Monte-Carlo 
method" as in Rogers's dual algorithm pifH]. 

In the next section, we briefly list (non-linear) Monte-Carlo algorithms which can be used to solve 
PDEs ([5])-([6| and highlight their weaknesses in the context of credit valuation adjustment. 

3. NON-LINEAR MONTE-CARLO ALGORITHMS 

3.1. A brute force algorithm. Using Feynman-Kac's formula, the solution of PDE ^ can be 

represented stochastically as 

(7) u{t,x) = e-^^^-*^EtA^{XT)]+ J (3e-^^'~*^Et,Au+is,Xs)]ds 

with X an Ito diffusion with generator £ and ]Et,a;[-] — E[-\Xt ~ x]. By assuming that the intensity 
/3 is small, we get the approximation (this is exact for PDE ([6]n ) 

(8) u(t,a;) = e-^(^-*)E,,,[V'(XT)]+/?e-^(^-*)y E,,,[(E,,xJ^(XT)])+]ds + 0(/?2) 
Then, at a next step, we discretise the Riemann integral 

n 

u{t,x) ~ e-^^^-*^Et,A^P{XT)]+Pe'^^^-*^Y.^tA{^u,x,M^T)])^]^U 

1=1 

This last expression can be numerically tackled by using a brute force "Monte-Carlo of Monte- 
Carlo" method. The second MC is used to compute Et.Xfi'i^T)] on each path generated by the 
first MC algorithm. Although straightforward, this method suffers from the curse of dimensionality 



^Precisely, we get e-^2(T-t)f:_^^^^^,^XT)] + Aa /f e-^2(»-t)Et,,[(l - R) (E,,x, [V-C^t)])"^ + HE,,x, [^{XtW-s. 



4 PIERRE HENRY-LABORDERE 

and requires generating 0{Ni x A'^2) paths. Due to this complexity, the hterature focuses on expo- 
sition of Hnear portfoHos for which the second MC can be skipped by using closed-form formulas 
or low-dimensional parametric regressions (see for example 3J in which the authors consider the 
pricing of CMS spread option and CCDSs). 

Could we design a simple (non-linear) Monte-Carlo algorithm which solves our PDEs ©-(IB]), 
without relying on an approximation such as (Is])? This is the purpose of this paper. 

3.2. Backward stochastic differential equations. A first approach is to simulate a backward 
stochastic differential equation (in short BSDE): 

(9) dXt = ^iit,Xt)dt + <j{t,Xt).dWu Xo=x 

(10) dYt = -PYt+dt + Zt(7(t,Xt).dWt 

(11) Yt = ^j{Xt) 

where {Y, Z) are required to be adapted processes and C = J^i Mi^x' + 5 J2i ji'^'^*)ij^%xi ■ BSDEs 



differ from (forward) SDKs in that we impose the terminal value (see Equation (111). Under the 
condition ip G L^{il), this BSDE admits a unique solution [TJ]. A straightforward application of 
Ito's lemma gives that the solution of this BSDE is {Yt = e^(^~*)u(t, Xt), Zt = e'^^'^-*'>Vxu{t, Xt)) 
with u the solution of PDE ^. This leads to a Monte-Carlo like numerical solution of ^ via an 
efficient discretization scheme for the above BSDE. 

This BSDE can be discretized by an Eulcr-like scheme {Yt-_^ is forced to be J"t;_^-adapted, {J^t)t>o 
being the natural filtration generated by the Brownian motions): 

Et,_jytJ - Yt^_, = -/3AU {eY+_^ + (1 - e)Et,_,[Yu]+) 

with 6 e [0, 1]. This is equivalent to (we take 6(3Ati < 1) 

This requires the computation of the conditional expectation Ef._jy(.] (in practise by regression 
methods) which could be quite difficult and time-consuming, especially for multi-asset portfolios. 

3.3. Gradient representation. A more powerful approach is synthesized by the following propo- 
sition which relies on Kunita's stochastic flows of diffeomorphisms (see [E]). Let u be the solution 
of the one- dimensional semi- linear PDE 

(12) dtu+la^it,x)dlu + fiu) = 



with the terminal condition u{T,x) = '4>{x). By differentiating equation (12 1 with respect to x 
(assuming smoothness of the coefficients) we get 

(13) 9tA + Uad^a) d, + ^a\t, x)dl) A + /' (u) A = 

with the terminal condition A{T,x) = 4''{x). The equation satisfied by the gradient A is then 
interpreted as a (linear) Fokker-Planck PDE. We have the following representation [16] 

u{t,x) = - f i)'{a)da Et[l(X^ - 2;)e-''t^ ■'^'("(^+*-"'-^°»'^''] 



where the Ito process X" is the solution to 

dX'^ = a{T + t- s, X^)dB, + {ad^a) (T + t-s, Xf)ds , s e [t,T] , X^ = a 
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Bg is a standard Brownian. This representation leads to a particle algorithm [16] . Although 
appealing, this (forward) approach is only applicable in the one-dimension setup for which we can 
use a PDE solver. Can we design a similar forward algorithm applicable in higher dimensions? 
This leads us to branching diffusions. 



3.4. Branching diffusions: an introduction. Branching diffusions have been first introduced 
by McKean jS^ to give a probabilistic representation of the Kolmogorov-Petrovskii-Piskunov PDE 
and more generally of semi-linear PDEs of the type 



(14) dtu + Cu + f3{t)(j2 Pku'' - u I = in 



X W^ 



\k=0 

u{T, x) = "0(2^) in K 



with /3(-) S M"*". Here the non-linearity is a power series in u where the coefficients satisfy the 
restrictive condition: 

C30 OO 

(15) f{u) = Y,Pku\ ^Pfc = l 0<pk<l 

k=a fc=o 

The probabilistic interpretation of such an equation goes as follows: Let a single particle start at 
the origin, perform an Ito diffusion on M.'^ with generator C, after a mean /?(•) exponential time 
(independent of X) die and produce k descendants with probability p^ {k = means that the 
particle dies without generating descendants). Then, the descendants perform independent Ito 
diffusions on E"* (with same generator C) from their birth locations, die and produce descendants 
after a mean /3(-) exponential times, etc. This process is called a d-dimcnsional branching diffusion 
with a branching rate /3(-). /3 can also depend spatially on x or be itself stochastic (Cox process). 

We note Zt = I z^^, . . . , z^ * j g M''^^' the locations of the particles alive at time t and Nt the 
number of particles at t (see Fig. IT for examples with 2 and 3 descendants). We consider then the 
multiplicative functional defined bj J 

Nt 



(16) u{t,x)=Et^,[l[i;{zir) 



j=i 



where Et,2;[-] = E[-|Aft — 1, z^ = x]. Note that as Nt can become infinite when m = X^feLo ^Pf^ -^ 
1 (super-critical regime, see }T2]). a sufficient condition on tp in order to have a well-behaved 



product is jV'l < 1- Then u solves the semi- linear PDE (14). This stochastic representation can be 



understood as follows: Mathematically, by conditioning on r, the first-time to jump of a Poisson 



process with intensity /3(t), we get from (161 



OO k Ni^(r) 

U{t, X) = Et,, [U>tV'(4)] + ^t,. [lr<T Y. Pf^^^ til n ^(-^ 

k=0 j = l i=l 



Yl ^ = 1 by convention. 
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where zl^'^^ is the position of the z-th particle at maturity T produced by the j-th particle generated 
at time r. By using the independence and the strong Markov property, we obtain 



\{t,x) = '^t,x\^r>T-^[^\)\ 



E 

fc=0 



iEt,.[ir<TPfcn^-[ n V'(4''"") 



k 



= Et,x[lr>TlPiz^)] + IEt^x[lr<T X^P*^ IT "(^' ""r)] 



k=0 



k=0 



E 



fc=0 



Then, by assuming that HV^Hoo < 1, u is uniformly bounded by 1 in [0,r] x M.'^ and we get from 
the Feynman-Kac formula that u is a viscosity solution to PDE ( 14) (see Theorem 6.4 in [TTj). By 
assuming that PDE (14) satisfies a comparison principle, we conclude that u = u. 



A first attempt in order to obtain a larger class of non-linearities than those defined by (151 is to 



consider an infinite collection of branching diffusions, the so-called super-diffusions. (15|) is then 
extended to 



(17) 



^{u) = au + bu^+ n{dr)[e 

Jo 



1 



where a > 0, 6 > and n is a Radon measure on (0, oo) satisfying L [r /\r'^)n(dr) < oo. The class 
of non-linearity as defined by (17) is more general than (15 1, in particular contains au -(- hv? with 



arbitrary positive coefficients a and h. Unfortunately, this requires a large number of branching 
diffusions (as the default intensity diverges) and the non-linearity is still restrictive. This leads us 
to introduce a new class of branching diffusions that can be traced back to Le Jan-Sznitman [9" in 
the context of stochastic (Fourier) representations of solutions of the incompressible Navier-Stokes 
equation. 



4. Marked branching diffusions 



The PDE ( 14 1 should be compared with the semi-linear PDE ^ arising in the pricing of coun- 
terparty risk. It seems too restrictive and unreasonable to approximate the non-linearity u^ by a 



polynomial of type ( 15 ) or even ( 17 1. A natural question is therefore to search if this construction 



can be generalized for an arbitrary polynomial for which the PDE is 
(18) dtu + Cu + P{F{u)-u)^Q 

with F[u) = X]fc=o ^kU^ an Af-order polynomial in u that we write for convenience F{u) — 
X]fe=o ( ~ ) Pku'^- We will show below that this can be achieved by counting the branching of each 
monomial u''. 



Assumption (Comp): In order to have uniqueness in the viscosity sense, we assume PDE (18) 
satisfies a comparison principle for sub- and super-solutions (see 0). 

For each Galton- Watson tree, we denote tuk G N the number of branching of monomial type u'^ 
with k € {0, . . . , M}. The descendants are drawn with an arbitrary distributi on p k - for example 
we can take a uniform distribution p^ — jfTj (see an other choice in section 4.3). In Fig. 1 we 
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LU 
Q 

Cl 



LU 
Q 

Cl 



Figure 1. Marked Galton- Watson random tree for the non-linearity F{u) 



P2' 



P2U + ^p^u'^ . The red (resp. blue) vertex corresponds to the weight ^ (resp. 
-). The diagram with two red vertices has the weights (wi = 2,a;2 = 0). 



have drawn the diagrams for the non-linearity F{u) 
then define the multiplicative functional: 



-P2M 



-p^u^ up to two defaults. We 



Main formula: 

(19) u(t,a;)=Et,,[nV'(^T)n 



Nt 



Pk 



Wfc = ttbranching type k 



4=1 k=0 

We state our main result (the proof is reported in the appendix): 

Theorem 4.1. Let us assume that u £ L°°([0,T] x M.'^) and ('Compj holds. The function u{t,x) 



is the unique viscosity solution of (18) 



Diagrammar interpretation. From Feynman-Kac's formula, we have 

(20) U{t, X) = Et,,[U>T^(XT)] + Et,,[F(u(r, Xr))lr<T] 

This integral equation can be recursively solved in terms of multiple exponential random times r^: 

u{t,x) = ¥.tAK,>THXT)] 

(21) + Et,,[i^(E,Jl,,>TV^(XT)] +E,„[F(E,Jl,,>T^(XT)])l.,<T])lro<T] + ••• 

Each term can be interpreted as a Feynman diagram (see Fig. [I]) representing the trajectory of a 
branching diffusion with a weight depending on the branching of each monomial. For example in 
Fig. [Tl the diagram with two red vertices corresponds to 



EtA^ro<T^rA^r^>Ti'{XTWro%,<T^r.A^r,>Ti^{XT)f 
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By assuming that the series (21 1 is convergent, one can guess that the solution is given by our 
multiphcative functional (19 1. 

In the next section, we focus on convergence issues and deduce a sufficient condition to ensure that 

u e L°°([0, T] X M'^) if V' is bounded. 

4.1. Convergence issues. The number of particles N{uj), produced by the branching uj = {ojq, . . . , ujm), 
is 

M 

(22) iVH = ^(fc-lK + l 

fc=0 

The probability of such a configuration satisfies the recurrence equation 

M T 

P(T|w)=/3^ / dtP{t\iUo,...,uJk~l,---,uJM)N{ujn,...,iJk-l,...,iJM) 

(23) p^g-/3fe(T-t)g-/3(T-t)(Ar(t^O:...,'^fc-l....,"M)-l) 

Indeed, if we have a tree with a branching {luq, . . . , Wfe — 1, . . . ,ujm) at time t, a particle among 
the N{uJo, . . . ,u)k — 1,---,wm) particles must die and produce k descendants (with probability 
Pklie.~^^^'^~^'). The remaining N{luq, . . . , a;^ — 1, . . . , ojm) ~ 1 particles must survive until maturity 
T (with probability e-/3(T'-t)(^('^o,...,t^fc-i,...,<^M)-i)). 

We prove in the appendix that the Laplace transform of P, P(T, c) — E[nfc=o e"''*"''], satisfies the 
equation 

/-P(T,c) , M 

M 

(25) P(T,c) = l if ^pfee-^'==l 

fc=0 

In the particular case of one branching type fc y^ 1, we have 



P(r,cfc) 



(l - e/^^C^-i) + e=''+'3'^('=^i)) '^ 



By assuming that -tp E L°°(M'^), the expectation in (19) can then be bounded by 



— ) IIV'II^^"^] = IIV'lloo'P f T -1. ^ -In ll,/,l|fe-l 
Pk J 

from which we deduce a sufficient condition for convergence 



(26) \u{0,x)\ < Eo,.[n — ||V'||^("^] = ||V'l|ooP T,-ln^-ln 



Proposition 4.2. Let us assume that ip e L°°(M''). Set p{s) = l3 {-s + J2tLo\ak\\\^\\^^s 



exists X G M^ such that 



(1) Case Y.k=o lafelllV'll^^ > 1- We have u G L°°([0,T] x W^) (as defined by ^^) if there 

'^ ds 

=T 

1 P{s) 

In the particular case of one branching type k, the sufficient condition for convergence reads 
as 



KIIHir(i-e-^^^'=-^0 



< 1 
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(2) Case Y.IU lafclllV-IlL"' < l-' « e L-([0,r] x W') for all T. 
Note that our blow-up criteria does not depend on the probabihties pk as expected. 

4.2. PDE Q. We assume that the function (1 — R)x'^ + Rx can be well approximated by a 
polynomial F{x) (see section [s]) and we consider the PDE 

dtu{t, x) + Cu{t, x) + -^ {FiEtAi^iXT)]) - u{t, x)) = 0, u{T, x) = ^j{x) 
From Feynman-Kac's formula, we have 

with T a Poisson default time with intensity /?/(! — R). As compared to the previous section, 
we have the term 'KtAF{Kr['4^{XT)])^T<T\ instead of Et^a;[i^(M(T, X^))1^<t]. This term can be 
computed using the previous algorithm by imposing that the particle can default only once. This 
corresponds to the first three diagrams in Fig. (II]). As Nt is valued in [0, M], our formula (19) is 
convergent here for all polynomial non-linearities. 

As a conclusion, without any modification, the branching particle algorithm can solve the two 
PDEs ((5])-((6| modulo that the non-linearly u^ can be fairly well approximated by a polynomial. 

4.3. Optimal probabilities pk- Is there a better choice than an uniform distribution pk — jfri 
for improving the convergence? 

For the PDE ([S]), the variance of the algorithm (depending on the probabilities pk) is bounded by 



(see Equation (26)) 



T,-21n^-21n||^||^-i 
Pk 



By minimizing with respect to p^ , we get 

Wk 



(27) Pk - M , 

Similarly, for the PDE (pi), the variance (depending on the probabilities pk) is bounded by 

M 2 



By minimizing with respect to pk, we get also (27) 



We recall that the population in the Galton- Watson tree disappears in finite time almost surely if 
m = X]/c=o ^Pk ^ 1 (see |12j). In the super-critical case m > 1, the population explodes at a finite 



time Texp with probability 1 — sq where sq = inf{s e [0, ^,J2k=oPkS = s}. From (27), we are in 
the super-critical case if J2k=o(^ ~ l)kfelll'/'llL > 0. 

4.4. Numerical Experiments. Before applying our algorithm to the problem of credit valuation 



adjustment, we check it on polynomials which do not belong to the classes defined by (15) and 

01. 



10 



PIERRE HENRY-LABORDERE 



N 


Fair(PDE2) 


Stdev(PDE2) 


Fair(PDEl) 


Stdev(PDEl) 


12 


20.78 


0.78 


21.31 


0.79 


14 


22.25 


0.39 


21.37 


0.39 


16 


21.97 


0.19 


21.76 


0.20 


18 


21.90 


0.10 


21.51 


0.10 


20 


21.86 


0.05 


21.48 


0.05 


22 


21.81 


0.02 


21.50 


0.02 



Table 1 . MC price quoted in percent as a function of the number of MC paths 
2^. PDE pricer(PDEl) = 21.82. PDE pricer(PDE2) = 21.50. Non-Unearity 
F{u)^Uu^-u^). 



N 


Fair(PDE2) 


Stdcv(PDE2) 


Fair(PDEl) 


Stdev(PDEl) 


12 


21.14 


0.78 


20.00 


0.78 


14 


21.56 


0.38 


19.90 


0.39 


16 


21.62 


0.19 


20.25 


0.20 


18 


21.31 


0.10 


20.39 


0.10 


20 


21.38 


0.05 


20.36 


0.05 


22 


21.36 


0.02 


20.40 


0.02 



Table 2. MC price quoted in percent as a function of the number of MC paths 
2^. PDE pricer(PDEl) = 21.37. PDE pricer(PDE2) = 20.39. Non-Unearity 



4.4.1. Experiment 1. We have implemented our algorithm for the two PDE types 

dtu + Cu + l3{F{u)-u) = Q, -u(r,2;) == l^>i : PDE2 
and 



dtu + Cu + p{F{¥.tA^XT>i]) - w) = 0, u{T, x) = U 



>i 



PDEl 



with F{u) ~ i {v?' — u^). C is the Ito generator of a geometric Brownian motion with a volatility 
cbs = 0.2 and the Poisson intensity is /3 = 0.05. In financial terms, this corresponds to a CDS 



spread around 500 basis points. The maturity is T = 10 years. From (27), we note that our optimal 



probability distributions for PDEl and PDE2 coincide with the uniform distribution. Moreover 



Proposition (4.2 1 gives that the solution does not blow up. 



The numerical method has been checked against a one-dimensional PDE solver with a fully implicit 
scheme (see Table, [l]) for which we find u = 21.82% (PDEl) and u = 21.50% (PDE2). Note 
that this algorithm converges as expected and the error is properly indicated by the Monte-Carlo 
standard deviation estimator (see column Stdev). 



4.4.2. Experiment 2. Same test with F{u) — ^ 
as above. 



*) (see Table. p| and same comments 



4.4.3. Experiment 3: Blow-up explosion. It is well-known that the semi-linear PDE in 



dtu + Cu 







blows up in finite time if d < 2 for any bounded positive payoff (see [H]). We deduce that the 
PDE with the non-linearity F{u) = v? +u blows up in finite time (Tmax) in one dimension. Using 
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Maturity(Year) 


BBM alg.(Stdev) 


PDE 


0.5 


71.66(0.09) 


71.50 


1 


157.35(0.49) 


157.17 


1.1 


oo(cx)) 


oo 



Table 3. MC price quoted in percent as a function of the maturity for the non- 
Hnearity F{u) ~ u^ + u. ipix) = la;>i. 



Proposition (4.2), our sufficient condition reads as 



We have verified this explosion when the maturity T is greater than 1 year (in our case iJj = la;>Oi 
1 1 "01 loo — 1) using our algorithm (and a PDE solver as a benchmark). Note that for T = 1, 
the algorithm starts to blow up (see Stdev = 0.49). A different stochastic representation can be 
obtained by setting u ~ e^-^^*^w. We get 

dtv + Cv + e^^-*)^^ -v = , v{T, x) = ^/j{x) 



and this can be interpreted as a binary tree with a weight e'"^ '^' 
gives then 



Our stochastic representation 



(28) 



Nt 

u{t, x) = e^-*E,^, [ [] 7/-(z^)e^"=r""" 



\T--r,) 



where Ti is the time where the i-ih. branching appears. This representation (28) appears in [TT] 
and was used to reproduce Sugitani's blow-up criteria |15) . 



5. Credit valuation adjustment algorithm 



In the previous section, we have assumed that the payoff was bounded: ?A € L°°. Then, the solution 

where v satisfies 



u can then be written as v 
(29) 



IIV-IU 



< 1 



dtv + Cv + 13 {v+ ~ v) ^ {), \\v{T, 

Therefore, by re-scaling, we can consider that the payoff satisfies the condition ||'0||oo < 1- The 
condition tp & L°° can be easily relaxed as observed in (|5], see Remark 3.7). Let V' be a payoff 
with a-exponential growth for some a > 0. We scale the solution by an arbitrary smooth positive 
function p given by 

p{x) = e"l^l for \x\ > M 
v{t,x) = p^^{x)v{t,x) 

If we write the linear operator C as Cv — fi{t,x)dxV + ^a^{t,x)d'^v, then v satisfies a PDEnwith 
the same non-linearity f3v^: 

dtv + Cv + I3{i+ ~d) ^0 
with Cv = (^p + cr^p^^dxp) dxV + ^cr^{t, x)d1v + [pp^^dxP + ^p^^a^dlp) v. 



What remains to be done in order to use (19) is to approximate v^ by a polynomial F{v) 
(30) dtV + Cv + P {F{v) - v) = 0, v{T,x)^ip{x) 



C is written in d = 1. A similar expression can be obtained in a multi-dimensional setup. 
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Figure 2. u^ versus its polynomial approximation on [—1. 1]. 



In our numerical experiments, we take (see Fig. ^ 

(31) F{u) = 0.0589 + 0.5u + 0.8164^^ - 0.4043m^ 



Proposition 4.2 gives that the solution does not blow up if /3T < 0.50829 (Take X = oo with 
IIV'lloo = !)• Moreover, as a numerical check of (26l, we have computed using a PDF solver the 
solution of ^ with il){x) = 1, F{u) = 0.0589 + 0.5w + 0.8164^2 + 0.4043u'', /3 = 0.05 and T = 10 



years. The solution X 



(32) 



T. - In i^ 
Pk 



X 



coincides with our upper bound in ( 26 ) and should satisfy 

ds 

0.5 



0.0589 + 0.5m + 0.8164^2 + 0.4043^4 



We found X = 4.497 (PDF solver) and the reader can check that this value satisfies the above 
identity (32) as expected. 



5.1. Algorithm: Final recipe. The algorithm for solving PDFs (l5])-(l6]) can be described by the 
following steps: 



^M 



(1) Choose a polynomial approximation of m+ ~ X]fc=o '^fcW on the domain [—1, 1]. 

(2) Simulate the assets and the Poisson default time with intensity j3 (resp. y^r) for PDF2 
(resp. PDFl). Note that the intensity /? can be stochastic (Cox process), usually calibrated 
to default probabilities implied from CDS market quotes. 



(3) At each default time, produce k descendants with probability pk (given by (27)). For PDF 
type 2, descendants, produced after the first default, become immortal. 
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(4) Evaluate for each particle alive the payoff 

n ^(4)( °^^ 'pf +^ ) n("^^^) (here,^.., = Oorl),PDEl 

where lo-^ denotes the number of branching type fc. We should highlight that the algo- 
rithm for PDEl is always convergent for all T whatever condition on the payoff as the 
multiplicative functional involves at most M particles. 

Remark 5.1. In the case of collateralized positions, the non-linearity uf should be substituted 
with {ut — Mt+A)"*^ where A is a delay. Using our polynomial approximation, we get F{ut — Mj+a). 
By expanding this function, we get monomials of the form {u^u^,^}. Our algorithm can then be 
easily extended to handle this case. At each default time r, we produce p descendants starting at 
(r, Xt) and q descendants starting at (r + A, Xt+a)- 

A natural question is to characterize the error of the algorithm as a function of the approximation 
error of u"*" by F{u). Using the parabolicity of the semi-linear PDE, we can characterize the bias 
of our algorithm (the proof is reported in the appendix): 

Proposition 5.2. Let us assume that F(v) and F{v) are two polynomials satisfying ('Compj, the 

sufficient condition in Prop. [JTEfor a maturity T and 

F{x) <x+< F{x) 



We denote v and v the corresponding solutions of (30) and v the solution of (29). Then 

V < V < V 



A similar result can be found for PDE ^. In the case of American options, our algorithm gives 
robust lower and upper bounds. 

5.2. Complexity. By approximating m"*" with an infinite high-order polynomial - say N2 - our 
algorithm converges towards the brute force "Monte-Carlo of Monte-Carlo" method with a com- 



plexity 0{Ni X N2). By comparison, with our choice (31|, the complexity is at most 0(4A^i) for 
PDE type (l6|. Moreover, this method allows to solve exactly PDE type ^, which can not be 
tackled without relying on an approximation within the "Monte-Carlo of Monte-Carlo" method. 

5.3. Numerical examples. We have implemented our algorithm for the two PDE types 



dtu + -x^ol^dlu + p(u+ ~u) ^ 0, u(T, a;) = 1 - 2.1^>i : PDEl 

and 

dtU + \x^alsdlu + ^^ ((1 - R)¥.tA^ - 2.1x^>i]+ + R&tA^ - 2.1xr>i] - u) = 0, PDE2 

with Poisson intensities j3 — 1%, j3 ~ 3% and a recovery rate R = 0.4 (see Tab. |4) [5J l6J l7| . In 
financial term, this corresponds to CDS spreads around 100 and 300 basis points. The method 



has been checked using a PDE solver with the polynomial approximation (31) (see Column "PDE 



with poly"). In order to justify the validity of (31 ), we have included the PDE price with the true 



non-linearity u"*" (see Column "PDE"). As it can be observed, prices, produced by our algorithm. 
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Maturity(Year) 


PDE with poly. 


BBM alg. 


PDE 


2 


11.62 


11.63(0.00) 


11.62 


4 


16.54 


16.53(0.00) 


16.55 


6 


20.28 


20.27(0.00) 


20.30 


8 


23.39 


23.38(0.00) 


23.41 


10 


26.11 


26.09(0.00) 


26.14 



Table 4. MC price quoted in percent as a function of the maturity for PDE 1 
with 13 = 1%. 



Maturity(Year) 


PDE with poly. 


BBM alg.(Stdev) 


PDE 


2 


11.62 


11.64(0.00) 


11.63 


4 


16.56 


16.55(0.02) 


16.57 


6 


20.32 


20.30(0.00) 


20.34 


8 


23.45 


23.45(0.00) 


23.48 


10 


26.20 


26.18(0.00) 


26.24 



Table 5. MC price quoted in percent as a function of the maturity for PDE 2 

with 13 = 1%. 



Maturity(Year) 


PDE with poly. 


BBM alg. 


PDE 


2 


12.34 


12.35(0.00) 


12.35 


4 


17.72 


17.71(0.00) 


17.75 


6 


21.77 


21.76(0.00) 


21.82 


8 


25.07 


25.06(0.00) 


25.14 


10 


27.89 


27.88(0.00) 


27.98 



Table 6. MC price quoted in percent as a function of the maturity for PDE 1 
with 13^3%. 



Maturity(Ycar) 


PDE with poly. 


BBM alg.(Stdcv) 


PDE 


2 


12.38 


12.39(0.00) 


12.39 


4 


17.88 


17.86(0.00) 


17.91 


6 


22.08 


22.07(0.01) 


22.14 


8 


25.58 


25.57(0.01) 


25.66 


10 


28.62 


28.60(0.01) 


28.74 



Table 7. MC price quoted in percent as a function of the maturity for PDE 2 
with 13 = 3%. 



converge to the PDE solver with the polynomial approximation and are close to the exact CVA 
values. We would like to highlight that replacing the Black-Scholes generator |a;^cr|g9^ by a 
multi-dimensional operator C can be easily handled in our framework by simulating the branching 
particles with a diffusion process associated to C. This is out-of-reach with finite-difference scheme 
methods and not such an easy step for the BSDE approach. 
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6. Conclusion 

Credit valuation adjustment is now an important quantitative issue which needs to receive special 
attention. The brute force "Monte-Carlo of Monte-Carlo" or the BSDE approach is not, as it looks 
like, a decent solution for multi-asset portfolios. We have shown the efficiency of our algorithm 
based on marked branching diffusions on various numerical examples. This method can also be 
used for semi-linear PDEs with polynomial non-linearities and extended to fully non-linear PDEs 
by including in the branching process Malliavin weights for derivatives. We left this investigation 
for future research. 

Acknowledgements. The author wishes to thank the members of the Global Markets Quantita- 
tive Research Group at Societe Generale for their comments. He is also grateful to Jean-Frangois 
Delmas and Denis Talay for useful discussions. 

Appendix 



Proof of Theorem 4-1 The proof proceeds similarly as in subsection 3.4 By using the indepen 



dence and the strong Markov property, we obtain 

M k N'^{t) m , -.uil 

U{t,x) = Et,4lr>Ti^{z^)]+Y.^t,^i^r<TakY[EA'[[ Y[[ — ] H4'''^)] 

fe=0 j=l i=l k=0 ^P'^^ 

M k 

= Et^a:[^r>Ttp{zT)] + Et,j;[lr<T X! "'^ IT ^^"^ ' ^r)] 

k=0 j = l 

= E,,,[e-/"^W'*^V'(4)]+/ Et^A/3is)e-^^'^^^'^''^F{uis,zl))]ds 

By assuming that u € L°°([0,T] x M"^), we deduce that w is a viscosity solution of PDE (18) 
(see Theorem 6.4 in \XL)- The comparison result (Assumption (Comp)) implies uniqueness, i.e. 

u = ii. D 

Proof of formula\24\ We set P(T|a;) = e~''"^^^"-'(7(T|w) for convenience. We get the relation 

M T 

q{T\uj) =PY. / ^^^9(^1^0, . . . , Wfc - 1, . . . , WAf )^(c^o, . . . , Wfc - 1, . . . , WMKe'^*('=-i^ 

which is equivalent to 

dTq{T\u) ^ /3^g(T|u;o, . . . ,c^fc - 1, . . . , c^a/)7V(c^o, . . • ,Wfc - 1, . . . , WM)Pfee'^^('=-i) 

fe=0 

q{Q\uj) = (5^=0 
The Laplace transform of q, q{T,c) = E«[nflo e^^''^^^'^"""], satisfies the first-order PDE 

M / M \ 

dTq{T\c) = /3^pJ q{T,c) - ^5,,g(r,c) e(/'^-^'=)('=-i) , g(0|c) = 1 

fc=0 \ 9=0 / 

The solution is given by 
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where the coefficients {cg(r)}^=o,....A/ are solutions of the ODEs 

fe=0 

The solution is given by Cq{t) — Cq ~ pt — \nU {t\c) with 

"^ = P i-Um + |:p.e(^^--)(^-i)C/'=(i|c)') , urn = 1 

This gives 

M 



q{T\c) = e^^U{T\c) if ^pfce('5^-^^)('=-i) ^ 1 

fe=0 
M 



/c=0 



where C/(T|c) satisfies 



s + EkLoPke^^^-'-^^"''^ 



1) c/e 



fc=0 



Finally, we use that P(r|c) = U{T\ ^ + /ST). D 

Proof of Proposition \5.S\ The function 6 = v — v satisfies the linear PDE 

dtS + C6 - pS + f3 { — l^-ij + (i (F(v) ~ v+) = , S(T, x)=0 

\ V ~ V J ^ ' 

Note that the term r^ = 1 — ( "j ~^' J Ivt^vt is lower bounded. Feynman-Kac's formula gives 

Sit, x)^ J (3EtA{Fiv) - v+) e-^ ^^' '^"'^"] 
from which we conclude the proof as F{x) > x^ by assumption. D 
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